namespace hms::cubic::norm_noLin
{

/* crude, but it's needed for both functions in exactly the same way,
 * and using functions would give overly complex return types,
 * because a std::array allows putting the 1XsMaps, but not the bufMap.
 * And the evaluation doesn't really belong into the same function */
#define CUBIC_ROOT_MAP_AND_EVAL \
	static_assert( std::is_base_of_v<Eigen::DenseBase<Buf>, Buf> ); \
	static_assert( std::is_base_of_v<Eigen::DenseBase<B>, B> );     \
	Index size { bXpr.size() };                                     \
	assert( calcBuf.size() >= bufCols * size );                     \
	ArrayXNsMap<bufCols> bufMap { calcBuf.data(), size, bufCols };  \
	Array1XsMap                                                     \
		b        { bufMap.col(0).data(), 1, size },                 \
		d        { bufMap.col(1).data(), 1, size },                 \
		discr    { bufMap.col(2).data(), 1, size },                 \
		p        { bufMap.col(3).data(), 1, size },                 \
		q        { bufMap.col(4).data(), 1, size };                 \
	b = bXpr;                                                       \
	d = dXpr;                                                       \
	p = cubic::norm_noLin::p(b  );                                  \
	q = cubic::norm_noLin::q(b,d);                                  \
	discr = cubic::discriminant(p,q);


template<typename B, typename D, typename K, typename Buf>
auto root(
	const B& bXpr, const D& dXpr,
	const K& k, 
	Buf& calcBuf
) -> decltype(auto)
{
	CUBIC_ROOT_MAP_AND_EVAL

	/* we didn't need this case yet, therefore YAGNI */
	if ( (discr > 0).any() ){
		throw std::runtime_error( "Unhandled discriminant" );
	}

	/* Buffer columns 1-4 can be reused at different points */
	Array1XsMap
		scale    { bufMap.col(1).data(), 1, size },
		acosTerm { bufMap.col(2).data(), 1, size },
		t        { bufMap.col(3).data(), 1, size };

	/* helper vars that are needed for all three roots */
	scale = cubic::threeReal::scale(p);
	acosTerm = cubic::threeReal::acosTerm(p,q);
	/* k decides which root to output */
	t = cubic::threeReal::depressedRoot( scale, acosTerm, k );
	return cubic::norm_noLin::root( t, b );
}

template<typename B, typename D, typename Buf>
auto root(
	const B& bXpr, const D& dXpr,
	Buf& calcBuf
) -> decltype(auto)
{
	CUBIC_ROOT_MAP_AND_EVAL

	/* make sure there is only a single real root */
	if ( (discr <= 0).any() ){
		throw std::runtime_error( "Unhandled discriminant" );
	}

	/* Buffer columns 1-4 can be reused at different points */
	Array1XsMap
		sqrtTerm { bufMap.col(1).data(), 1, size },
		t        { bufMap.col(2).data(), 1, size };

	sqrtTerm = cubic::singleReal::sqrtTerm(p,q);
	t        = cubic::singleReal::depressedRoot(q, sqrtTerm);

	return cubic::norm_noLin::root( t, b );
}

#undef CUBIC_ROOT_MAP_AND_EVAL

} // namespace hms::cubic::norm_noLin
